Riquesa total

riquesa_total<- dades_finals %>% 
  summarise(riquesa_total=n_distinct(ID))
riquesa_total

## # A tibble: 1 × 1
##   riquesa_total
##           <int>
## 1           137

Riquesa per cada mostreig i tractament

riquesa <- dades_finals %>%
  group_by(Mostreig, Tractament) %>%
  summarise(riquesa = n_distinct(ID))

## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.

riquesa

## # A tibble: 11 × 3
## # Groups:   Mostreig [6]
##    Mostreig             Tractament riquesa
##    <chr>                <chr>        <int>
##  1 Campus-Juliol-2021   NS              20
##  2 Campus-Juliol-2021   S               21
##  3 Campus-Jun-2021      NS              45
##  4 Campus-Jun-2021      S               34
##  5 Campus-juliol-2020   NS              39
##  6 Campus-juliol-2020   S               24
##  7 Campus-maig-2020     NS              45
##  8 Franqueses-Maig-2021 NS              20
##  9 Franqueses-Maig-2021 S               11
## 10 Sabadell-Juny-2021   NS              16
## 11 Sabadell-Juny-2021   S               12

Riquesa per cada zona i tractament

riquesa_zones <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(riquesa = n_distinct(ID))

## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.

riquesa_zones

## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament riquesa
##   <chr>      <chr>        <int>
## 1 Campus     NS              98
## 2 Campus     S               59
## 3 Franqueses NS              20
## 4 Franqueses S               11
## 5 Sabadell   NS              16
## 6 Sabadell   S               12

Riquesa per grups

# Abelles
riquesa_abelles<- dades_finals %>% 
  filter(Ordre == 'HymenopteraAbella') %>% 
  summarise(riquesa_abelles=n_distinct(ID))
riquesa_abelles

## # A tibble: 1 × 1
##   riquesa_abelles
##             <int>
## 1              80

# Coleòpters
riquesa_coleoptera<- dades_finals %>% 
  filter(Ordre == 'Coleoptera') %>% 
  summarise(riquesa_coleoptera=n_distinct(ID))
riquesa_coleoptera

## # A tibble: 1 × 1
##   riquesa_coleoptera
##                <int>
## 1                 33

# Sírfids
riquesa_sirfids<- dades_finals %>% 
  filter(Familia == 'Syrphidae') %>% 
  summarise(riquesa_sirfids=n_distinct(ID))
riquesa_sirfids

## # A tibble: 1 × 1
##   riquesa_sirfids
##             <int>
## 1               9

Abundància total

abundancia_total<-dades_finals %>% 
  summarise(Numero_total=n())
abundancia_total

## # A tibble: 1 × 1
##   Numero_total
##          <int>
## 1         1111

Abundància per cada zona i tractament

abundancia_zona <- dades_finals %>%
  group_by(Zona, Tractament) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.

abundancia_zona

## # A tibble: 6 × 3
## # Groups:   Zona [3]
##   Zona       Tractament Numero_total
##   <chr>      <chr>             <int>
## 1 Campus     NS                  649
## 2 Campus     S                   323
## 3 Franqueses NS                   45
## 4 Franqueses S                    14
## 5 Sabadell   NS                   29
## 6 Sabadell   S                    51

Abundància de cada gènere total

ab_genere <- dades_finals %>%
  group_by(Genere) %>%
  summarise(Numero_total = n())
head(ab_genere)

## # A tibble: 6 × 2
##   Genere     Numero_total
##   <chr>             <int>
## 1 Acmaeodera            6
## 2 Amegilla              1
## 3 Andrena               6
## 4 Anthaxia             29
## 5 Anthidium            10
## 6 Apis                 11

Abundància de cada ID per zona i tractament

abundaciazona_ID <- dades_finals %>%
  group_by(Zona, Tractament, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Zona', 'Tractament'. You can override
## using the `.groups` argument.

head(abundaciazona_ID)

## # A tibble: 6 × 4
## # Groups:   Zona, Tractament [1]
##   Zona   Tractament ID                     Numero_total
##   <chr>  <chr>      <chr>                         <int>
## 1 Campus NS         Acmaeodera cylindrica             5
## 2 Campus NS         Amegilla sp.                      1
## 3 Campus NS         Anthaxia (Anthaxia)               1
## 4 Campus NS         Anthaxia (Melanthaxia)            5
## 5 Campus NS         Anthaxia hungarica               10
## 6 Campus NS         Anthaxia sp.                      7

Abundancia de cada ID per cada mostreig i tractament

abundacia_ID <- dades_finals %>%
  group_by(Mostreig, Tractament, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.

head(abundacia_ID)

## # A tibble: 6 × 4
## # Groups:   Mostreig, Tractament [1]
##   Mostreig           Tractament ID                   Numero_total
##   <chr>              <chr>      <chr>                       <int>
## 1 Campus-Juliol-2021 NS         Amegilla sp.                    1
## 2 Campus-Juliol-2021 NS         Anthaxia sp.                    1
## 3 Campus-Juliol-2021 NS         Bombus pascuorum                1
## 4 Campus-Juliol-2021 NS         Ceratina cucurbitina            1
## 5 Campus-Juliol-2021 NS         Ceratina parvula                7
## 6 Campus-Juliol-2021 NS         Chrysididae                     1

Corba d’acumulació d’espècies per mostreig (11 mostrejos)

corba <- specaccum(data_shannon_MT, permutations = 500)

## Warning in cor(x > 0): the standard deviation is zero

plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MT)

##     Species     chao  chao.se jack1 jack1.se    jack2     boot  boot.se  n
## All     137 236.8148 31.40658   207 26.07681 250.1182 167.3704 12.88879 11

plot(poolaccum(data_shannon_MT))

En aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 236 (chao), 207 (jackknife primer ordre), 250 (jackknife segon ordre), i 167 (bootstrap).

Corba d’acumulació d’espècies segons mostreig, tractament i trampa

corba <- specaccum(data_shannon_MTTR, permutations = 500)
plot(corba)

Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)

specpool(data_shannon_MTTR)

##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se   n
## All     137 224.1493 29.51259 201.3564 11.27322 241.7798 164.5183 5.987207 101

plot(poolaccum(data_shannon_MTTR))

Amb aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 224 (chao), 201 (jackknife primer ordre), 241 (jackknife segon ordre), i 164 (bootstrap).

Corba i extrapolació pel Campus:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades pel Campus
abundacia_trampa_campus <- dades_finals %>%
  filter(Zona == 'Campus') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.

#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxcampus <- cast(abundacia_trampa_campus, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxcampus <- as.data.frame(taula_auxcampus)
taula_auxcampus[is.na(taula_auxcampus)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_campus <- taula_auxcampus %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_campus) <- data_shannon_campus$Agrupacio
data_shannon_campus <- data_shannon_campus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivcampus <- diversity(data_shannon_campus)
print(shannondivcampus)

##  Campus-juliol-2020_NS_2  Campus-juliol-2020_NS_3  Campus-juliol-2020_NS_4 
##                1.7917595                1.8310205                2.0449312 
##  Campus-juliol-2020_NS_6  Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10 
##                1.7328680                1.8310205                1.0735428 
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13 
##                0.0000000                1.7478681                1.6434177 
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16 
##                1.1240357                1.5498260                1.3862944 
## Campus-juliol-2020_NS_20   Campus-juliol-2020_S_1   Campus-juliol-2020_S_5 
##                2.2718685                0.6931472                1.9072840 
##   Campus-juliol-2020_S_7   Campus-juliol-2020_S_9  Campus-juliol-2020_S_17 
##                0.6924212                0.0000000                1.0986123 
##  Campus-juliol-2020_S_18  Campus-juliol-2020_S_19  Campus-juliol-2020_S_21 
##                1.6094379                1.5607104                1.1189689 
##  Campus-juliol-2020_S_22  Campus-Juliol-2021_NS_2  Campus-Juliol-2021_NS_4 
##                1.0751393                0.0000000                1.3862944 
##  Campus-Juliol-2021_NS_6  Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10 
##                0.9502705                1.4648164                0.6931472 
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16 
##                1.3296613                1.3321790                0.0000000 
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20   Campus-Juliol-2021_S_1 
##                1.3862944                1.0986123                1.0986123 
##   Campus-Juliol-2021_S_3   Campus-Juliol-2021_S_5   Campus-Juliol-2021_S_7 
##                0.6365142                2.0963526                1.5000764 
##   Campus-Juliol-2021_S_9  Campus-Juliol-2021_S_11  Campus-Juliol-2021_S_13 
##                1.1973401                1.4388830                0.0000000 
##  Campus-Juliol-2021_S_15  Campus-Juliol-2021_S_17     Campus-Jun-2021_NS_2 
##                0.6931472                0.5091373                2.3345491 
##     Campus-Jun-2021_NS_4     Campus-Jun-2021_NS_6     Campus-Jun-2021_NS_8 
##                2.0854684                1.1537419                2.6197294 
##    Campus-Jun-2021_NS_10    Campus-Jun-2021_NS_12    Campus-Jun-2021_NS_14 
##                1.8310205                2.5575077                1.7623137 
##    Campus-Jun-2021_NS_16    Campus-Jun-2021_NS_18    Campus-Jun-2021_NS_20 
##                1.3593685                1.9061547                0.9502705 
##      Campus-Jun-2021_S_1      Campus-Jun-2021_S_3      Campus-Jun-2021_S_5 
##                2.2718685                2.0028830                1.7351265 
##      Campus-Jun-2021_S_7      Campus-Jun-2021_S_9     Campus-Jun-2021_S_11 
##                0.9502705                2.1383972                1.5595812 
##     Campus-Jun-2021_S_13     Campus-Jun-2021_S_15     Campus-Jun-2021_S_17 
##                1.8576598                0.7356219                2.1639557 
##     Campus-Jun-2021_S_19    Campus-maig-2020_NS_1    Campus-maig-2020_NS_2 
##                1.9722470                0.8392696                1.3642812 
##    Campus-maig-2020_NS_3    Campus-maig-2020_NS_4    Campus-maig-2020_NS_5 
##                1.6957425                2.5521719                1.5171064 
##    Campus-maig-2020_NS_6    Campus-maig-2020_NS_7    Campus-maig-2020_NS_8 
##                2.2614951                1.2424533                1.7480673 
##    Campus-maig-2020_NS_9   Campus-maig-2020_NS_10   Campus-maig-2020_NS_11 
##                1.6769878                1.9722470                1.7480673 
##   Campus-maig-2020_NS_12   Campus-maig-2020_NS_13   Campus-maig-2020_NS_14 
##                1.0986123                1.5481155                1.8891592 
##   Campus-maig-2020_NS_15   Campus-maig-2020_NS_16   Campus-maig-2020_NS_17 
##                1.8392967                1.0356097                1.3296613 
##   Campus-maig-2020_NS_18   Campus-maig-2020_NS_19   Campus-maig-2020_NS_20 
##                1.5810938                1.3535914                0.9002561

corbacampus <- specaccum(data_shannon_campus, permutations = 500)
plot(corbacampus)

# Extrapolació de les dades
specpool(data_shannon_campus)

##     Species     chao  chao.se    jack1 jack1.se    jack2     boot  boot.se  n
## All     118 184.0553 24.55884 170.3457 9.341038 201.8116 140.7347 5.098803 81

Corba i extrapolació per Franqueses:

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Franqueses
abundacia_trampa_franqueses <- dades_finals %>%
  filter(Zona == 'Franqueses') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.

#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxfranq <- cast(abundacia_trampa_franqueses, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxfranq <- as.data.frame(taula_auxfranq)
taula_auxfranq[is.na(taula_auxfranq)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_franqueses <- taula_auxfranq %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_franqueses) <- data_shannon_franqueses$Agrupacio
data_shannon_franqueses <- data_shannon_franqueses %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivfranq <- diversity(data_shannon_franqueses)
print(shannondivfranq)

## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3 
##                 1.0789922                 0.8675632                 0.6931472 
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5  Franqueses-Maig-2021_S_1 
##                 2.2739657                 0.7549968                 1.3321790 
##  Franqueses-Maig-2021_S_2  Franqueses-Maig-2021_S_4 
##                 1.0986123                 1.5607104

corbafranqueses <- specaccum(data_shannon_franqueses, permutations = 500)
plot(corbafranqueses)

# Extrapolació de les dades
specpool(data_shannon_franqueses)

##     Species     chao  chao.se  jack1 jack1.se    jack2     boot  boot.se n
## All      26 52.32292 16.96108 42.625 8.605049 53.01786 33.13315 4.330973 8

Corba i extrapolació per Sabadell

#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Sabadell
abundacia_trampa_sabadell <- dades_finals %>%
  filter(Zona == 'Sabadell') %>%
  group_by(Mostreig, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.

#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxsaba <- cast(abundacia_trampa_sabadell, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxsaba <- as.data.frame(taula_auxsaba)
taula_auxsaba[is.na(taula_auxsaba)] <- 0

# Agrupar Mostreig amb Tractament i Trampa
data_shannon_sabadell <- taula_auxsaba %>%
  mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
  ungroup() %>%
  dplyr::select(-Mostreig, -Tractament, -Trampa)

# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_sabadell) <- data_shannon_sabadell$Agrupacio
data_shannon_sabadell <- data_shannon_sabadell %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)

# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivsaba <- diversity(data_shannon_sabadell)
print(shannondivsaba)

##  Sabadell-Juny-2021_NS_1  Sabadell-Juny-2021_NS_4  Sabadell-Juny-2021_NS_7 
##                1.7478681                0.9502705                1.7917595 
##  Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11   Sabadell-Juny-2021_S_2 
##                1.4750763                1.3862944                0.6365142 
##   Sabadell-Juny-2021_S_3   Sabadell-Juny-2021_S_5   Sabadell-Juny-2021_S_6 
##                0.5623351                1.4995094                1.0986123 
##   Sabadell-Juny-2021_S_9  Sabadell-Juny-2021_S_10  Sabadell-Juny-2021_S_12 
##                0.6931472                0.5004024                1.2852930

corbasabadell <- specaccum(data_shannon_sabadell, permutations = 500)
plot(corbasabadell)

# Extrapolació de les dades
specpool(data_shannon_sabadell)

##     Species     chao  chao.se    jack1 jack1.se   jack2     boot  boot.se  n
## All      23 45.45833 17.10731 35.83333 5.316379 44.4697 28.44918 2.631169 12

T-student per comparar la diversitat de Shannon entre mostrejos i tractaments

#Convertir la taula de shannondiv en dataframe
shannon_ttest <- as.data.frame(shannondiv)

# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
shannon_ttest$row <- rownames(shannon_ttest)

# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
shannon_ttest <- shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)

# Prova t-student i boxplots
shannonT <- shannon_ttest %>% 
  group_by(Mostreig) %>% 
  summarise(
    mean_S = mean(shannondiv[Tractament=='S']),
    mean_NS = mean(shannondiv[Tractament=='NS']),
    t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(shannondiv ~ Tractament)$p.value,
                             NA),.groups = "drop")
  
shannonT

## # A tibble: 6 × 4
##   Mostreig              mean_S mean_NS t_test_result
##   <chr>                  <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1.02    0.964        0.844 
## 2 Campus-Jun-2021        1.74    1.86         0.637 
## 3 Campus-juliol-2020     1.08    1.54         0.0840
## 4 Campus-maig-2020     NaN       1.56        NA     
## 5 Franqueses-Maig-2021   1.33    1.13         0.566 
## 6 Sabadell-Juny-2021     0.897   1.47         0.0232

ggplot(shannon_ttest, aes(x = factor(Mostreig), y = shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon x Zona i Tractaments",
     x = "Mostreig",
     y = "Diversitat Shannon")

Només en el mostreig de Sabadell la diversitat de Shannon és significativament diferent (p valor < 0.05): l’índex de Shannon és major al tractament NS.

Corba d’acumulació. Riquesa en funció del nº de trampes (specaccum)

# A la taula data_shannon_MTTR extreure una nova columna que indiqui el mostreig
 metadata_shannon <- data_shannon_MTTR %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) 

# View(metadata_shannon)

# Subset each mostreig into its own dataframe (i eliminar la columna de "mostreig" perquè totes les dades de la taula siguin numèriques)
metadata_shannon %>% filter(Zona == 'Campus') %>% dplyr::select(-Zona)-> Campus
metadata_shannon %>% filter(Zona == 'Franqueses') %>% dplyr::select(-Zona) -> Franqueses
metadata_shannon %>% filter(Zona == 'Llig') %>% dplyr::select(-Zona) -> Llig
metadata_shannon %>% filter(Zona == 'Mura') %>% dplyr::select(-Zona) -> Mura
metadata_shannon %>% filter(Zona == 'Sabadell') %>% dplyr::select(-Zona) -> Sabadell
metadata_shannon %>% filter(Zona == 'TV') %>% dplyr::select(-Zona) -> TV7041

# Corba d'acumulació d'espècies per cada mostreig
corba1 = specaccum(Campus, method = "random")
corba2 = specaccum(Franqueses, method = "random")
corba3 = specaccum(Sabadell, method = "random")

# Creating data frames for each set of data
data1 <- data.frame(Sites = corba1$sites, Richness = corba1$richness, SD = corba1$sd)
data1$Zona <- "Campus"

data2 <- data.frame(Sites = corba2$sites, Richness = corba2$richness, SD = corba2$sd)
data2$Zona <- "Franqueses"

data3 <- data.frame(Sites = corba3$sites, Richness = corba3$richness, SD = corba3$sd)
data3$Zona <- "Sabadell"


# Combine all data frames into one
dades_combinades <- rbind(data1, data2, data3)

library(ggplot2)

# Gráfico de dispersión con líneas y error bars
ggplot(dades_combinades, aes(x = Sites, y = Richness, color = Zona)) +
  geom_point() +
  geom_line(aes(group = Zona), alpha = 0.7) +  # Líneas
  # geom_ribbon(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD, fill = Mostreig), alpha = 0.01) +  # Incertidumbre
  geom_errorbar(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD), width = 0.2, alpha = 0.5) +  # Barras de error
  ylim(0,NA)+
  labs(x = "Nº trampes", y = "Riquesa", color = "Zona") +
  theme_minimal()

Corba de rarefracció. Riquesa en funció del nombre d’individus (rarefy)

dades_agrupades <- metadata_shannon %>%
  group_by(Zona) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades) 

# Buscar el nombre mínim d'observacions en una zona (és 84, a TV juliol)
abundacia_trampa <- dades_finals %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.

abundacia_trampa

## # A tibble: 564 × 7
## # Groups:   Zona, Mes, Any, Tractament, Trampa [101]
##    Zona   Mes      Any Tractament Trampa ID                         Numero_total
##    <chr>  <chr>  <dbl> <chr>       <dbl> <chr>                             <int>
##  1 Campus Juliol  2021 NS              2 Lasioglossum glabriusculum            1
##  2 Campus Juliol  2021 NS              4 Curculionidae                         1
##  3 Campus Juliol  2021 NS              4 Lasioglossum malachurum               1
##  4 Campus Juliol  2021 NS              4 Pompilidae                            1
##  5 Campus Juliol  2021 NS              4 Pyronia cecilia                       1
##  6 Campus Juliol  2021 NS              6 Lasioglossum glabriusculum            1
##  7 Campus Juliol  2021 NS              6 Lasioglossum politum                  3
##  8 Campus Juliol  2021 NS              6 Pompilidae                            1
##  9 Campus Juliol  2021 NS              8 Chrysididae                           1
## 10 Campus Juliol  2021 NS              8 Hylaeus sp.                           1
## # ℹ 554 more rows

min_n <- abundacia_trampa %>% 
  group_by(Zona) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades) <- dades_agrupades$Zona
dades_agrupades <- dades_agrupades[,-1]
dades_agrupades <- dades_agrupades %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades,sample=min_n)

##     Campus Franqueses   Sabadell 
##   28.12144   26.00000   19.50913 
## attr(,"Subsample")
## [1] 59

rarecurve(dades_agrupades,step=5)

# Mirar vídeo i posar-ho més maco?

Corba d’interpolació: riquesa en funció del nombre d’individus (iNEXT). (??)

library(iNEXT)
transposades <- t(dades_agrupades)
D_abund <- iNEXT (transposades, datatype = 'abundance')
plot (D_abund)

Corbes de rarefracció d’abelles per cada tractament

# Extreure la columna de tractament a la taula d'abundàncies per ID de cada mostreig
dades_shannon_abelles <- data_shannon_abelles %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])

# Agrupar per tractament
dades_agrupades_tractament <- dades_shannon_abelles %>%
  group_by(Tractament) %>%
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

#View (dades_agrupades_tractament) 

dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abelles <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.

#View (abundacia_trampa_abelles)

min_n <- abundacia_trampa_abelles %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dades_agrupades_tractament) <- dades_agrupades_tractament$Tractament
dades_agrupades_tractament <- dades_agrupades_tractament[,-1]
dades_agrupades_tractament <- dades_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dades_agrupades_tractament,sample=min_n)

##       NS        S 
## 59.89088 37.00000 
## attr(,"Subsample")
## [1] 321

rarecurve(dades_agrupades_tractament,step=5)

Corbes de rarefracció d’abelles al Campus per tractament

# Extreure la columna Zona de dades_shannon_abelles i filtrar per Campus
dades_shannon_abellescampus <- dades_shannon_abelles %>% 
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
dades_shannon_abellescampus <- dades_shannon_abellescampus[dades_shannon_abellescampus$Zona == "Campus", ]

# Agrupar per tractament
dadescampus_agrupades_tractament <- dades_shannon_abellescampus %>%
  group_by(Tractament) %>%
  summarise(across(where(is.numeric), sum)) %>% 
  as.data.frame()

#View (dadescampus_agrupades_tractament) 

# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abellescampus <- dades_abelles %>%
  group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
  summarise(Numero_total = n())

## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.

#View (abundacia_trampa_abellescampus)

min_n <- abundacia_trampa_abellescampus %>% 
  group_by(Tractament) %>% 
  summarize (n_especimens = sum(Numero_total)) %>% 
  summarize(min=min(n_especimens)) %>% 
  pull(min)

rownames(dadescampus_agrupades_tractament) <- dadescampus_agrupades_tractament$Tractament
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament[,-1]
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament %>% 
                    mutate_all(as.numeric)

vegan::rarefy(dadescampus_agrupades_tractament,sample=min_n)

## Warning in vegan::rarefy(dadescampus_agrupades_tractament, sample = min_n):
## requested 'sample' was larger than smallest site maximum (266)

##       NS        S 
## 54.01949 35.00000 
## attr(,"Subsample")
## [1] 321

rarecurve(dadescampus_agrupades_tractament,step=5)

Riquesa total

riquesa_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     4.33    3.1         0.303 
## 2 Campus-Jun-2021        7.3     9.5         0.190 
## 3 Campus-juliol-2020     3.89    5.92        0.0281
## 4 Campus-maig-2020     NaN       6.55       NA     
## 5 Franqueses-Maig-2021   4       4.6         0.743 
## 6 Sabadell-Juny-2021     3.43    4.8         0.201

ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Generalment es compleix un patró de major riquesa al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0,0281)

Riquesa d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  

## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
riquesa_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_abelles <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Riquesa ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abelles

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     3.89    2.1        0.0847 
## 2 Campus-Jun-2021        4.3     5.4        0.338  
## 3 Campus-juliol-2020     2.78    4.46       0.00961
## 4 Campus-maig-2020     NaN       3.84      NA      
## 5 Franqueses-Maig-2021   2       3          0.341  
## 6 Sabadell-Juny-2021     3.14    3.8        0.568

 ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Riquesa d'abelles per Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Generalment es compleix un patró de major riquesa d’abelles al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0.0096)

Riquesa de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  

## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
riquesa_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 

  proves_t_coleoptera <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_coleoptera

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        2.88    3.7          0.357
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       2.44        NA    
## 5 Franqueses-Maig-2021   1.5     1            0.500
## 6 Sabadell-Juny-2021     1       1.33         0.423

  ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Riquesa de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  

## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
riquesa_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Riquesa = n_distinct(ID),.groups = "drop") 
  
  proves_t_sirfids <- riquesa_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Riquesa[Tractament == "S"]),
      mean_NS = mean(Riquesa[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Riquesa ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_sirfids

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA

ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Riquesa x Zona i Tractaments",
       x = "Mostreig",
       y = "Riquesa")

Hi ha poques dades per obtenir resultats del t-test.

Abundància total

abundancia_x_trampa <- dades_finals %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abtotal <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_abtotal

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    14.1     4.5          0.130
## 2 Campus-Jun-2021       11.6    18.5          0.113
## 3 Campus-juliol-2020     8.89    8.54         0.904
## 4 Campus-maig-2020     NaN      15.4         NA    
## 5 Franqueses-Maig-2021   4.67    9            0.178
## 6 Sabadell-Juny-2021     7.29    5.8          0.606

 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

No hi ha cap resultat significatiu (p-valor<0.05). No hi ha un patró clar de com varia l’abundància en funció del tractament.

Abundància d’abelles

# Filtrar les dades d'abelles
  dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles)  

## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
abundancia_x_trampa <- dades_abelles %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_ababelles <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = ifelse(length(unique(Tractament)) == 2, 
                             t.test(Abundancia ~ Tractament)$p.value,
                             NA)
      ,.groups = "drop")
  
  proves_t_ababelles

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021    13.7     3.5         0.108 
## 2 Campus-Jun-2021        7.6     8.9         0.589 
## 3 Campus-juliol-2020     7.44    7.08        0.896 
## 4 Campus-maig-2020     NaN       8.68       NA     
## 5 Franqueses-Maig-2021   2       4.4         0.0902
## 6 Sabadell-Juny-2021     7       4.8         0.454

 ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

A Franqueses l’abundància d’abelles és significativament superior al tractament NS (p-valor=0.09015). En la resta de mostrejos el patró no és clar.

Abundància de coleòpters

# Filtrar les dades per coleòpters
  dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera)  

## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1        29 Campus     Juliol         7  2021      4 NS         <NA>            
## 2       158 Campus     Juliol         7  2021     10 NS         <NA>            
## 3       183 Campus     Juliol         7  2021      5 S          <NA>            
## 4       101 Franqueses Maig           5  2021      2 NS         <NA>            
## 5       150 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## 6       151 Franqueses Maig           5  2021      5 NS         Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
abundancia_x_trampa <- dades_coleoptera %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n(),.groups = "drop") 
  
  proves_t_abcoleoptera <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_abcoleoptera

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl>         <dbl>
## 1 Campus-Juliol-2021     1       1           NA    
## 2 Campus-Jun-2021        4.12    9.2          0.175
## 3 Campus-juliol-2020     1.33    1.25         0.851
## 4 Campus-maig-2020     NaN       6.61        NA    
## 5 Franqueses-Maig-2021   2       4.67         0.496
## 6 Sabadell-Juny-2021     1       1.33         0.423

  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.

Abundància de sírfids

# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids)  

## # A tibble: 6 × 26
##   Especimen Zona       Mes    `Mes num`   Any Trampa Tractament Identificacio   
##       <dbl> <chr>      <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>           
## 1       226 Campus     Juliol         7  2021      7 S          <NA>            
## 2       265 Campus     Juliol         7  2021     14 NS         <NA>            
## 3       277 Franqueses Maig           5  2021      4 NS         <NA>            
## 4       278 Franqueses Maig           5  2021      4 NS         <NA>            
## 5       279 Franqueses Maig           5  2021      4 NS         <NA>            
## 6       280 Franqueses Maig           5  2021      4 NS         Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Proves t-student i plots
abundancia_x_trampa <- dades_sirfids %>%
    group_by(Mostreig, Tractament, Trampa) %>%
    summarise(Abundancia = n_distinct(ID),.groups = "drop") 
  
  proves_t_absirfids <- abundancia_x_trampa %>%
    group_by(Mostreig) %>%
    summarise(
      mean_S = mean(Abundancia[Tractament == "S"]),
      mean_NS = mean(Abundancia[Tractament == "NS"]),
      t_test_result = tryCatch(
      t.test(Abundancia ~ Tractament)$p.value,
      error = function(e) NA
    )
      ,.groups = "drop")
  
  proves_t_absirfids

## # A tibble: 6 × 4
##   Mostreig             mean_S mean_NS t_test_result
##   <chr>                 <dbl>   <dbl> <lgl>        
## 1 Campus-Juliol-2021        1    1    NA           
## 2 Campus-Jun-2021           1    1    NA           
## 3 Campus-juliol-2020        3    1    NA           
## 4 Campus-maig-2020        NaN    1.18 NA           
## 5 Franqueses-Maig-2021      1    3    NA           
## 6 Sabadell-Juny-2021      NaN    1    NA

  ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
  geom_boxplot(width = 0.7, outlier.shape = NA) +
  facet_wrap(~Mostreig, scales = "free") +
  labs(title = "Boxplot Abundancia x Zona i Tractaments",
       x = "Mostreig",
       y = "Abundancia")

No hi ha prous dades per obtenir bons resultats de la prova t-student.

Anàlisi de similitud (ANOSIM) Total

ano = anosim(NMDS_data, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
ano

## 
## Call:
## anosim(x = NMDS_data, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.06667 
##       Significance: 0.721 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor > 0.05): la composició de les comunitats no és significativament diferent entre els tractaments S i NS. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.

ADONIS total

# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata_TOTAL <- NMDS_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)

Zona_Tractament_metadata_TOTAL

##                               Zona Tractament
## Campus-juliol-2020_NS       Campus         NS
## Campus-juliol-2020_S        Campus          S
## Campus-Juliol-2021_NS       Campus         NS
## Campus-Juliol-2021_S        Campus          S
## Campus-Jun-2021_NS          Campus         NS
## Campus-Jun-2021_S           Campus          S
## Campus-maig-2020_NS         Campus         NS
## Franqueses-Maig-2021_NS Franqueses         NS
## Franqueses-Maig-2021_S  Franqueses          S
## Sabadell-Juny-2021_NS     Sabadell         NS
## Sabadell-Juny-2021_S      Sabadell          S

NMDStot.dist <- vegdist(wisconsin(sqrt(NMDS_data)), k=2)

adonis(formula = NMDStot.dist~Tractament/Zona, data = Zona_Tractament_metadata_TOTAL, permutations = 1000)$aov.tab

## 'adonis' will be deprecated: use 'adonis2' instead

## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                 Df SumsOfSqs MeanSqs F.Model      R2  Pr(>F)  
## Tractament       1    0.3611 0.36111  1.1373 0.10191 0.25874  
## Tractament:Zona  4    1.5947 0.39868  1.2556 0.45004 0.03996 *
## Residuals        5    1.5877 0.31753         0.44805          
## Total           10    3.5435                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha diferències significatives entre tractaments (p-valor = 0.25874), però la interacció entre zona i tractament explica un 45% de la variabilitat, i és un resultat significatiu (p-valor=0.03996).

Espècies indicadores dels tractaments S i NS total

indicadores = multipatt(NMDS_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives de cap dels dos tractaments.

summary(indicadores, alpha = 0.2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 137
##  Selected number of species: 6 
##  Number of species associated to 1 group: 6 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  3 
##                    stat p.value
## Anthaxia sp.      0.422   0.182
## Eucera elongatula 0.422   0.181
## Ceratina parvula  0.391   0.178
## 
##  Group S  #sps.  3 
##                             stat p.value
## Vespula germanica          0.500   0.186
## Lasioglossum glabriusculum 0.430   0.184
## Lasioglossum griseolum     0.385   0.184
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Anthaxia sp., Eucera elongatula, Ceratina parvula), i 3 espècies indicadores de S (Vespula germanica, Lasioglossum gabriusculum i Lasioglossum griseolum). Però realment no son significatives.

Anàlisi de similitud (ANOSIM) Campus

ano = anosim(NMDS_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)

## 'nperm' >= set of all permutations: complete enumeration.

## Set of permutations < 'minperm'. Generating entire set.

ano

## 
## Call:
## anosim(x = NMDS_data_campus, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.03704 
##       Significance: 0.42857 
## 
## Permutation: free
## Number of permutations: 5039

No signiticatiu (p valor = 0.42857): la composició de les comunitats de S i NS al Campus no és significativament diferent entre els dos tractaments (S i NS).

ADONIS Campus

# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_CAMPUS <- NMDS_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

Tractament_metadata_CAMPUS

##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS

NMDScamp.dist <- vegdist(wisconsin(sqrt(NMDS_data_campus)), k=2)

adonis(formula = NMDScamp.dist~Tractament, data = Tractament_metadata_CAMPUS, permutations = 1000)$aov.tab

## 'adonis' will be deprecated: use 'adonis2' instead

## Set of permutations < 'minperm'. Generating entire set.

## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.27747 0.27747 0.87044 0.14828 0.5924
## Residuals   5   1.59387 0.31877         0.85172       
## Total       6   1.87134                 1.00000

El tractament només explica un 14,8% de la variabilitat, i el resultat no és significatiu.

Espècies indicadores dels tractaments S i NS Campus

indicadores = multipatt(NMDS_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 137
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives de cap dels dos tractaments.

summary(indicadores, alpha = 0.2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 137
##  Selected number of species: 5 
##  Number of species associated to 1 group: 5 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  3 
##                     stat p.value
## Halictus scabiosae 0.775   0.149
## Anthaxia sp.       0.542   0.147
## Ceratina parvula   0.498   0.147
## 
##  Group S  #sps.  2 
##                             stat p.value
## Lasioglossum glabriusculum 0.731   0.118
## Vespula germanica          0.707   0.145
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Halictus scabiosae, Anthaxia sp., i Ceratina parvula), i 2 espècies indicadores de S (Lasioglossum glabriusculum i Vespula germanica). Però no son resultats significatius.

Anàlisi de similitud (ANOSIM) abelles totals

ano = anosim(NMDSabelles_data, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
ano

## 
## Call:
## anosim(x = NMDSabelles_data, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.04267 
##       Significance: 0.6117 
## 
## Permutation: free
## Number of permutations: 9999

No signiticatiu (p valor > 0.6037): la composició de les comunitats és pràcticament idèntica (valor de l’estadístic R proper a 0, i fins i tot negatiu). A més, el resultat no és significatiu. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.

ADONIS abelles totals

# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata <- NMDSabelles_data %>%
  mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>% 
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)

Zona_Tractament_metadata

##                               Zona Tractament
## Campus-juliol-2020_NS       Campus         NS
## Campus-juliol-2020_S        Campus          S
## Campus-Juliol-2021_NS       Campus         NS
## Campus-Juliol-2021_S        Campus          S
## Campus-Jun-2021_NS          Campus         NS
## Campus-Jun-2021_S           Campus          S
## Campus-maig-2020_NS         Campus         NS
## Franqueses-Maig-2021_NS Franqueses         NS
## Franqueses-Maig-2021_S  Franqueses          S
## Sabadell-Juny-2021_NS     Sabadell         NS
## Sabadell-Juny-2021_S      Sabadell          S

NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data)), k=2)

adonis(formula = NMDSabe.dist~Tractament/Zona, data = Zona_Tractament_metadata, permutations = 1000)$aov.tab

## 'adonis' will be deprecated: use 'adonis2' instead

## Permutation: free
## Number of permutations: 1000
## 
## Terms added sequentially (first to last)
## 
##                 Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## Tractament       1    0.3748 0.37477  1.3134 0.11721 0.1299  
## Tractament:Zona  4    1.3958 0.34896  1.2229 0.43656 0.0979 .
## Residuals        5    1.4267 0.28534         0.44622         
## Total           10    3.1973                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El tractament només explica un 11,72% de la variabilitat observada entre les poblacions de S i NS, i no és un resultat significatiu. La interacció entre el tractament i la zona explica un 43,65% de la variabilitat observada, i és un resultat marginalment significatiu (0,0989, < 0.1).

Espècies d’abelles indicadores dels tractaments S i NS

indicadores = multipatt(NMDSabelles_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 80
##  Selected number of species: 0 
##  Number of species associated to 1 group: 0 
## 
##  List of species associated to each combination: 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives.

summary(indicadores, alpha = 0.2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 3 
##  Number of species associated to 1 group: 3 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                    stat p.value
## Eucera elongatula 0.422   0.181
## 
##  Group S  #sps.  2 
##                             stat p.value
## Lasioglossum glabriusculum 0.430   0.184
## Lasioglossum griseolum     0.385   0.184
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si augmentem el nivell de significació (alfa) a 0.2, hi ha 1 espècie indicadora a NS (Eucera elongatula) i 2 espècies indicadores a S (L. glabriusculum i L. griseolum). Però no son resultats significatius.

Anàlisi de similitud (ANOSIM) abelles Campus

ano = anosim(NMDSabelles_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)

## 'nperm' >= set of all permutations: complete enumeration.

## Set of permutations < 'minperm'. Generating entire set.

ano

## 
## Call:
## anosim(x = NMDSabelles_data_campus, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.03704 
##       Significance: 0.34286 
## 
## Permutation: free
## Number of permutations: 5039

La composició de les comunitats entre els dos tractaments és molt similar (valor de l’estadístic d’R molt proper a 0), i el resultat no és significatiu.

ADONIS abelles Campus

# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS1_Tractament_metadata <- NMDSabelles_data_campus %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

ADONIS1_Tractament_metadata

##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS

NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data_campus)), k=2)

adonis(formula = NMDSabe.dist~Tractament, data = ADONIS1_Tractament_metadata, permutations = 1000)$aov.tab

## 'adonis' will be deprecated: use 'adonis2' instead

## Set of permutations < 'minperm'. Generating entire set.

## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.26326 0.26326  0.9175 0.15505 0.4905
## Residuals   5   1.43466 0.28693         0.84495       
## Total       6   1.69792                 1.00000

El factor Tractament explica un 15.5% de la variabilitat d’abelles al Campus entre els tractaments S i NS, però aquesta explicació no és significativa.

Espècies d’abelles indicadores dels tractaments S i NS al Campus

indicadores = multipatt(NMDSabelles_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 2 
##  Number of species associated to 1 group: 2 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                     stat p.value
## Halictus scabiosae 0.775   0.149
## 
##  Group S  #sps.  1 
##                             stat p.value
## Lasioglossum glabriusculum 0.731   0.118
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No hi ha espècies indicadores significatives. Si augmenta el nivell de significació alfa = 0.2, sí que apareixen Halictus scabiosae com a indicador del tractament NO segat, i Lasioglossum glabriusculum com a indicador del tractament segat.

Anàlisi de similitud ANOSIM presència-absència

ano = anosim(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)

## 'nperm' >= set of all permutations: complete enumeration.

## Set of permutations < 'minperm'. Generating entire set.

ano

## 
## Call:
## anosim(x = NMDS_data_abellescampus_presabs, grouping = Tractament_metadata$Tractament,      permutations = 9999, distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.1019 
##       Significance: 0.62857 
## 
## Permutation: free
## Number of permutations: 5039

Les composicions de les comunitats dels dos tractaments son molt similars, i el resultat no és significatiu.

ADONIS presència-absència

# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS_Tractament_metadata <- NMDS_data_abellescampus_presabs %>%
  mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)

ADONIS_Tractament_metadata

##                       Tractament
## Campus-juliol-2020_NS         NS
## Campus-juliol-2020_S           S
## Campus-Juliol-2021_NS         NS
## Campus-Juliol-2021_S           S
## Campus-Jun-2021_NS            NS
## Campus-Jun-2021_S              S
## Campus-maig-2020_NS           NS

NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDS_data_abellescampus_presabs)), k=2)

adonis(formula = NMDSabe.dist~Tractament, data = ADONIS_Tractament_metadata, permutations = 1000)$aov.tab

## 'adonis' will be deprecated: use 'adonis2' instead

## Set of permutations < 'minperm'. Generating entire set.

## Permutation: free
## Number of permutations: 5039
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Tractament  1   0.23205 0.23205 0.97825 0.16363 0.3716
## Residuals   5   1.18603 0.23721         0.83637       
## Total       6   1.41808                 1.00000

El tractament explica un 16% de la variabilitat entre les comunitats, però el resultat no és significatiu.

Espècies indicadores presència-absència

indicadores = multipatt(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.2
## 
##  Total number of species: 80
##  Selected number of species: 1 
##  Number of species associated to 1 group: 1 
## 
##  List of species associated to each combination: 
## 
##  Group NS  #sps.  1 
##                     stat p.value
## Halictus scabiosae 0.775   0.137
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Només Halictus scabiosae com a indicadora del tractament NS, però el resultat no és significatiu (nivell de significació alfa = 0.2)

Model lineal mixt (tractament com a variable explicativa, i zona com a variable aleatòria)

# install.packages("lme4")
# install.packages("Matrix")

model_lasioglossum <- lmer(Abundancia_Total ~ Tractament + (1|Zona), data = lasioglossums_zona_tractament)
summary(model_lasioglossum)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Abundancia_Total ~ Tractament + (1 | Zona)
##    Data: lasioglossums_zona_tractament
## 
## REML criterion at convergence: 642.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1440 -0.4778 -0.1408  0.1962  7.2804 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept)  1.557   1.248   
##  Residual             35.226   5.935   
## Number of obs: 101, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)    1.926      1.172  2.051   1.643  0.23908   
## TractamentS    3.954      1.226 98.906   3.224  0.00171 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.429

plot(model_lasioglossum)

Quan el tractament és NS, l’abundància de Lasioglossums és superior a 0 (Intercept: marginalment significatiu). Quan el tractament és S, s’espera que l’abundància de lasioglossums augmenti 3,954 respecte el tractament NS, i és un resultat molt significatiu (p-valor = 0.00171). Els Lasioglossums son més abundants a les zones segades (concorda amb algunes de les espècies “indicadores” de S!)

Model lineal (Tractament com a variable explicativa, sense considerar la zona)

# Convertir la variable Tractament a factor
lasioglossums_zona_tractament$Tractament <- factor(lasioglossums_zona_tractament$Tractament)

# Ajustar el model lineal
model <- lm(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
summary(model)

## 
## Call:
## lm(formula = Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.421 -2.619 -1.421  1.381 43.579 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.619      0.754   3.473 0.000764 ***
## TractamentS    3.802      1.229   3.093 0.002576 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.985 on 99 degrees of freedom
## Multiple R-squared:  0.08811,    Adjusted R-squared:  0.0789 
## F-statistic: 9.566 on 1 and 99 DF,  p-value: 0.002576

# Crear el gráfico
plot(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament, pch = 16)

Abelles totals, amb Zona com a factor aleatori

#Obrir la base de dades amb la columna de les mesures de la ITS i filtrar dades per abelles
abelles_ITS <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]

head(abelles_ITS)

## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Model lineal mixt ITS (mida)
abelles_ITS_filtrades <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]
model_abelles_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|Zona), data = abelles_ITS_filtrades)
summary(model_abelles_ITS)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona)
##    Data: abelles_ITS_filtrades
## 
## REML criterion at convergence: 1453.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3928 -0.4191 -0.2937  0.1489  5.4311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Zona     (Intercept) 0.03073  0.1753  
##  Residual             0.48378  0.6955  
## Number of obs: 683, groups:  Zona, 3
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.78796    0.11569   2.44978  15.455  0.00167 ** 
## TractamentS  -0.51521    0.05431 676.60529  -9.487  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.217

La ITS de les abelles disminueix -0.51521 quan el tractament és S (abelles més petites quan està segat). És significatiu (p-valor = 2e-16)

# Plot del model
# Prediccions del model sobre les dades:
predicted_values <- fitted(model_abelles_ITS)
observed_values <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]$`Inter_tegular_span`

# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values, Predicted = predicted_values)

# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") +  # Línea de referencia de 45 grados
  labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
  theme_minimal()

ggplot(abelles_ITS_filtrades, aes(Tractament, Inter_tegular_span, fill = Tractament)) +
  geom_boxplot() +
  facet_wrap(~Zona) +
  labs(x = "Tratament", y = "ITS") +
  scale_fill_manual(values = c("S" = "#ADD8E6", "NS" = "#FFC0CB"))

Abelles del Campus, amb data_mostreig com a factor aleatori

# Filtrar les dades només per abelles del Campus
abellescampus_ITS <- abelles_ITS[abelles_ITS$Zona == "Campus",]

#Definir nova variable: "data_mostreig", combinació de mes i any de mostreig
abellescampus_ITS <- abellescampus_ITS %>%
  mutate(data_mostreig = paste(Mes, Any, sep = "-"))

head(abellescampus_ITS)

## # A tibble: 6 × 27
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>

# Fer el model
model_abellescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = abellescampus_ITS)
summary(model_abellescampus_ITS)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
## REML criterion at convergence: 1196.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3101 -0.5834 -0.2132  0.1382  5.4645 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.04629  0.2151  
##  Residual                  0.42778  0.6540  
## Number of obs: 594, groups:  data_mostreig, 4
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.54523    0.11430   3.29968  13.519 0.000529 ***
## TractamentS  -0.36675    0.06534 537.59881  -5.613 3.19e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.243

Molt significatiu: la ITS de les abelles (mida) disminueix -0.36675 quan el tractament és S. Les abelles son més petites quan està segat.

# Fer el plot
# Prediccions del model sobre les dades:
predicted_values1 <- fitted(model_abellescampus_ITS)
observed_values1 <- abellescampus_ITS[!is.na(abellescampus_ITS$Inter_tegular_span),]$`Inter_tegular_span`

# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values1, Predicted = predicted_values1)

# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") +  
  labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
  theme_minimal()

### ATENCIÓ!!! NO ESTIC SEGURA DE SI AQUESTA REPRESENTACIÓ GRÀFICA DEL MODEL ÉS ACURADA. SON VALORS PREDITS I OBSERVATS, O SON ELS RESIDUS AJUSTATS I NOSEQ???? VEURE GRÀFICS DELS ALTRES MODELS ANTERIORS.

ITS abelles per espècies totals:

# Crear nova base de dades: només una espècie per cada tractament.
especies_ITS <- abelles_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
head(especies_ITS)

## # A tibble: 6 × 26
##   Especimen Zona   Mes    `Mes num`   Any Trampa Tractament Identificacio
##       <dbl> <chr>  <chr>      <dbl> <dbl>  <dbl> <chr>      <chr>        
## 1         6 Campus Juliol         7  2021     18 NS         <NA>         
## 2         8 Campus Juliol         7  2021     18 NS         <NA>         
## 3         9 Campus Juliol         7  2021     18 NS         <NA>         
## 4        11 Campus Juliol         7  2021     18 NS         <NA>         
## 5        27 Campus Juliol         7  2021      4 NS         <NA>         
## 6       159 Campus Juliol         7  2021     10 NS         <NA>         
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## #   `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## #   NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## #   `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## #   Nesting_type <chr>, Lecty <chr>, Mostreig <chr>

# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especies_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especies_ITS)
summary(model_especies_ITS)

## 
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especies_ITS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3679 -0.7225 -0.3465  0.6686  3.5533 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.1199     0.1603  13.225   <2e-16 ***
## TractamentS  -0.2894     0.2581  -1.121    0.266    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.051 on 68 degrees of freedom
##   (35 observations deleted due to missingness)
## Multiple R-squared:  0.01815,    Adjusted R-squared:  0.003711 
## F-statistic: 1.257 on 1 and 68 DF,  p-value: 0.2662

No significatiu

ITS abelles per espècies del campus

# Crear nova base de dades: només una espècie per cada tractament.
especiescampus_ITS <- abellescampus_ITS %>%
  distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
View(especiescampus_ITS)

# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especiescampus_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especiescampus_ITS)
summary(model_especiescampus_ITS)

## 
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especiescampus_ITS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2394 -0.6734 -0.3183  0.6081  2.8422 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9914     0.1667  11.949   <2e-16 ***
## TractamentS  -0.2301     0.2582  -0.891    0.377    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9859 on 58 degrees of freedom
##   (33 observations deleted due to missingness)
## Multiple R-squared:  0.01351,    Adjusted R-squared:  -0.003501 
## F-statistic: 0.7941 on 1 and 58 DF,  p-value: 0.3765

No significatiu.

Nest_type en funció del tractament de les abelles totals (zona = factor aleatori)

abelles_ITS<-abelles_ITS[!is.na(abelles_ITS$Nesting_type),]

abelles_ITS["Nius"] = 0
abelles_ITS$Nius[abelles_ITS["Nesting_type"]=="S"]="S"
abelles_ITS$Nius[abelles_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia1 <- table(abelles_ITS$Tractament, abelles_ITS$Nius)
print(taula_contingencia1)

##     
##      Not_S   S
##   NS    71 324
##   S     14 291

# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia1)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Substrats de nidificació per tractament")

# Model:
abelles_ITS$Nius <- factor(abelles_ITS$Nius)
# Nomes Zona o tambe mostreig, com a factors aleatoris
model_abelles_nest <- glmer(Nius ~ Tractament + (1|Zona) + (1|Mostreig), data = abelles_ITS, family = binomial)

## boundary (singular) fit: see help('isSingular')

summary(model_abelles_nest)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nius ~ Tractament + (1 | Zona) + (1 | Mostreig)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    493.3    511.5   -242.6    485.3      696 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7820  0.2091  0.2282  0.4715  0.5146 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Mostreig (Intercept) 3.795e-02 0.1947978
##  Zona     (Intercept) 1.277e-10 0.0000113
## Number of obs: 700, groups:  Mostreig, 6; Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4843     0.1645   9.024  < 2e-16 ***
## TractamentS   1.6262     0.3450   4.714 2.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.442
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Significatiu! Les abelles del tractament S tenen més probabilitat (1.6262) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 2.43e-06).

# Plot del model
plot(fitted(model_abelles_nest), resid(model_abelles_nest), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Nest_type en funció del tractament de les abelles del Campus (data_mostreig = factor aleatori)

abellescampus_ITS<-abellescampus_ITS[!is.na(abellescampus_ITS$Nesting_type),]

abellescampus_ITS["Nius"] = 0
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]=="S"]="S"
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]!="S"]="Not_S"

# Taula de contingència
taula_contingencia2 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Nius)
print(taula_contingencia2)

##     
##      Not_S   S
##   NS    63 291
##   S     13 243

# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia2)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
  ggtitle("Substrats de nidificació per tractament")

# Model:
abellescampus_ITS$Nesting_type <- factor(abellescampus_ITS$Nius)
model_abellescampus_nest <- glmer(Nesting_type ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_nest)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Nesting_type ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    440.2    453.4   -217.1    434.2      607 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5162  0.2291  0.2360  0.4689  0.4997 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 0.02712  0.1647  
## Number of obs: 610, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4968     0.1764   8.484  < 2e-16 ***
## TractamentS   1.5004     0.3860   3.887 0.000101 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.512

Significatiu! Les abelles del tractament S tenen més probabilitat (1.5004) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 0.000101).

# Plot del model
plot(fitted(model_abellescampus_nest), resid(model_abellescampus_nest), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lecty en funció del tractament de les abelles totals (Zona = factor aleatori)

# Taula de contingència
taula_contingencia3 <- table(abelles_ITS$Tractament, abelles_ITS$Lecty)
print(taula_contingencia3)

##     
##      Oligolectic Polylectic
##   NS         142        249
##   S           38        263

# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia3)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty per tractament")

# Model:
abelles_ITS$Lecty <- factor(abelles_ITS$Lecty)
model_abelles_lecty <- glmer(Lecty ~ Tractament + (1|Zona), data = abelles_ITS, family = binomial)

summary(model_abelles_lecty)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | Zona)
##    Data: abelles_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    730.0    743.6   -362.0    724.0      689 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0678 -0.6616  0.3260  0.7371  1.5115 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Zona   (Intercept) 0.7457   0.8635  
## Number of obs: 692, groups:  Zona, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.2956     0.5350   0.553    0.581    
## TractamentS   1.6319     0.2235   7.302 2.84e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.132

Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat molt significatiu (p-valor=), la qual cosa suggereix un fort efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abelles_lecty), resid(model_abelles_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")

Lecty en funció del tractament de les abelles del campus (data_mostreig = factor aleatori)

# Taula de contingència
taula_contingencia4 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Lecty)
print(taula_contingencia4)

##     
##      Oligolectic Polylectic
##   NS         127        223
##   S           20        232

# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia4)

ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(x = "Tractament", y = "Nº", fill = "Lecty") +
  ggtitle("Lecty per tractament")

# Model:
abellescampus_ITS$Lecty <- factor(abellescampus_ITS$Lecty)
model_abellescampus_lecty <- glmer(Lecty ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)

summary(model_abellescampus_lecty)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Lecty ~ Tractament + (1 | data_mostreig)
##    Data: abellescampus_ITS
## 
##      AIC      BIC   logLik deviance df.resid 
##    468.8    482.0   -231.4    462.8      599 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2618  0.0518  0.1706  0.5796  1.2245 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  data_mostreig (Intercept) 6.291    2.508   
## Number of obs: 602, groups:  data_mostreig, 4
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   2.2170     1.3182   1.682   0.0926 .
## TractamentS   0.6375     0.3271   1.949   0.0513 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TractamentS -0.091

Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat marginalment significatiu (p-valor = 0.0513), la qual cosa suggereix un cert efecte del tractament en la lecticitat de les abelles.

# Plot del model
plot(fitted(model_abellescampus_lecty), resid(model_abellescampus_lecty), 
     xlab = "Valors ajustats", ylab = "Residus estandarditzats")